home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / dsp / fft / fft_eyal.lha / fft_eyal / fftg.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-31  |  12.0 KB  |  362 lines
  1. /* ------------------------------ fftg.c ------------------------------------ */
  2. /*                                                                            */
  3. /* Author:      Eyal Lebedinsky                                               */
  4. /* Date:        May 1990                                                      */
  5. /* Version:     9 June 1991                                                   */
  6. /*                                                                            */
  7. /* Program generator for integer DFFT. This program was developed for a very  */
  8. /* specific need and no attempt was made to make it general. Still, the basic */
  9. /* algorithm is a standard one from the literature. For maximum performance   */
  10. /* of the integer calculation a special mechanism was employed to track shift */
  11. /* operations that could be delayed. Coefficients are kept as 16 bits with 1  */
  12. /* sign bit, 1 integer bit and 14 fractional bits. So a typicall multiply     */
  13. /* will do 16bit*16bits=32bits, then shift down one bit and use the HIGH 16   */
  14. /* bits as the result. Before adding two numbers they must be scaled down one */
  15. /* bit to avoid overflow. Adding more numbers calls for more scaling. The     */
  16. /* delayed scaling allows the intermediate array to hold numbers that are     */
  17. /* scaled to a varying degree until a clash occurs. Well, this space is too   */
  18. /* short to elaborate any further and if you did similar kind of work then    */
  19. /* you should know what I am talking about.                                   */
  20. /*                                                                            */
  21. /* usage: fftgc File Size EpName                                              */
  22. /*      File - output file name.                                              */
  23. /*      Size - sample order: 8 for 256 points, 10 for 1024...                 */
  24. /*      EpName - name of the generated function ('_' added)                   */
  25. /* Now assemble the output file and then use as:                              */
  26. /*                                                                            */
  27. /* #define M 8                  [or whatever size]                            */
  28. /* #define N (1 << M)                                                         */
  29. /* short x[N];                  [input/output array]                          */
  30. /* short qf[(N/2)+1];           [power spectrum]                              */
  31. /* extern void far fft ()                                                     */
  32. /*                                                                            */
  33. /* [then call it as: fft();]                                                  */
  34. /*                                                                            */
  35. /* This is a modified version of Sorensen et al, IEEE ASSP-35 no 6.           */
  36. /*                                                                            */
  37. /* - uses 3+3 complex mult.                                                   */
  38. /* - no pre scrambling. descrambling done when power spectrum calculated.     */
  39. /* - moved to 0-origin as proper for c.                                       */
  40. /* - trig values in table.                                                    */
  41. /* - the case for sin==cos treated separately.                                */
  42. /* - uses delayed scaling down.                                               */
  43. /* - does length 256 fft with 642 mults.                                      */
  44. /*                                                                            */
  45. /*                                                                            */
  46. /* This program is released into the public domain.                           */
  47. /*                                                                            */
  48. /*----------------------------------------------------------------------------*/
  49.  
  50. #include <stdio.h>
  51. #include <math.h>
  52.  
  53. #define    MAXPOINTS    1024
  54. #define SC15 32768                /* scaling factor 2^15 */
  55.  
  56. /***
  57. About SC15:
  58.     Fractions are kept as 16 bit quantities with one sign bit, no integer
  59.     bits and 15 fractional bits. SC15 is used to convert a floating point
  60.     number (in the range -1 .. 1) into such an integer.
  61.  
  62.     There is no integer part because all fractions are really kept already
  63.     scaled down by a factor of 2. This is a result of the fact that all
  64.     fractions are constants (sin() and cos() values) and whenever we
  65.     multiply a point by a fraction we process to scale the result down by
  66.     one bit anywat - so the reduce constants save us some processing time.
  67. ***/
  68.  
  69. static short mm[MAXPOINTS+1];        /* digit reverse indices */
  70.  
  71. /***
  72. About mm[]:
  73.     The values in the sample points array x[] are not scramled at any
  74.     stage. Instead, the array mm[] indicated where a point should be
  75.     if we DID do the scrambling. All indices into x[] use mm[] and
  76.     this is how the power spectrum extracts the data. This means that
  77.     you CANNOT use the x[] array after calling fft() because the points
  78.     are not in place. If you want to change the output of fft() then just
  79.     modify the last part of this program (the power spectrum logic) to
  80.     calculate whatever you need.
  81. ***/
  82.  
  83. static int cntx[9];        /* stats counters */
  84.  
  85. /***
  86. About cntx[]:
  87.     This array gathers statistics:
  88.         0 memory loads
  89.         1 memory stores
  90.         2 register shifts
  91.         3 adds
  92.         4 mults
  93.         5 memory shifts due to smad() correction
  94. ***/
  95.  
  96. static char *file_name, *ep_name;
  97.  
  98. extern void    start_fft (), end_fft ();
  99. extern void    fft1 (),                /* scaling */
  100.         fft2 (), fft3 (), fft4 (), fft5 (),    /* SRFFT */
  101.         fft7 (), fft8 ();            /* power spectrum */
  102.  
  103. /***
  104. about ad[], mad and smad():
  105.     Because we do integer arithmetic, adding two numbers can overflow the
  106.     integer magnitude. For this reason the integers are scaled down (by one
  107.     bit) before the addition. Sometimes we add more that two numbers and it
  108.     is necessary to scale further.
  109.  
  110.     The FFT algorithm does not use all of the points in the same way in one
  111.     pass, which means that some data needs to be scaled and some does not
  112.     need to. Actualy, in the first pass some points are not used at all.
  113.     This means that SOME of the points would have been scaled down. If we
  114.     add (in a later stage) such a point with an un-scaled point the result
  115.     is invalid. For this purpose the array ad[] keeps track of how far each
  116.     point was scaled, smad() checks that a point is properly scaled and mad
  117.     is the needed scaling at this stage. As the FFT progresses, some points
  118.     are kept at a different scaling untill they are used in an expression
  119.     with other points. If smad() finds a point out of adjustment the it
  120.     issues a scaling instruction using the fft1() function.
  121.  
  122.     The program was adjusted to optimize the scaling (the various
  123.     fft2() ... fft5() have internal scaling rules) and you may notice that
  124.     smad() hardly ever needs to do any correction (usually there is only one
  125.     final shift). This makes the whole ad[] almost redundant but it was left
  126.     in for future use.
  127.  
  128.     If your smaple points are less that 16 bits then you still need to keep
  129.     scaling you data. You need to have enough spare bits to avoid the
  130.     scaling, and the number depends on the size of the sample - the larger
  131.     the sample the more stages the FFT has and the more spare bits you
  132.     need.
  133. ***/
  134.  
  135. static short ad[MAXPOINTS+1];    /* how far each point is adjusted */
  136. static int mad;         /* maximum adjust at this level */
  137.  
  138. static void
  139. smad (i, j)                    /* set proper adjust           */
  140. int    i, j;
  141. {
  142.     while (ad[i] < (mad-j)) {
  143.         fprintf (stderr, "x[%u] >>= %d (%d-%d)\n", i, 1, mad, j);
  144.         fft1 (mm[i]);
  145.         ++ad[i];
  146.         cntx[0] += 1;
  147.         cntx[1] += 1;
  148.         cntx[2] += 1;
  149.         cntx[5] += 1;
  150.     }
  151. }
  152.  
  153. main (argc, argv)
  154. int    argc;
  155. char    *argv[];
  156. {
  157.     int n, m, j, i, k, is, id, n1, n2, n4, n8, ind,
  158.         i0, i1, i2, i3, i4, i5, i6, i7, i8,
  159.         cc1, ss1, cc3, ss3;
  160.     float e, a, a3;
  161.  
  162. /* ---------- Handle arguments ------------------------------------------- */
  163.  
  164.     if (argc < 4) {
  165.         printf ("usage: fftg File PowerOf2 EpName [stderrFile]\n");
  166.         exit (0);
  167.     }
  168.  
  169.     file_name = argv[1];
  170.     sscanf (argv[2], "%u", &m);
  171.     ep_name = argv[3];
  172.     if (argc > 4)
  173.         freopen (argv[4], "wt", stderr);
  174.  
  175. /* ---------- Initialise  ------------------------------------------------ */
  176.  
  177.     start_fft (file_name, ep_name);
  178.  
  179.     n = 1 << m;
  180.  
  181.     for (i = 0; i < 9; ++i)
  182.         cntx[i] = 0;
  183.  
  184.     mad = 0;
  185.     for (i = 1; i <= n; ++i)
  186.         ad[i] = 0;
  187.  
  188. /* ---------- Digit reverse counter -------------------------------------- */
  189.  
  190.     for (i = 1; i <= n; ++i) {
  191.         mm[i] = i-1;        /* '-1' for shifting to 0-origin    */
  192.     }
  193.  
  194.     j = 1;
  195.     n1 = n - 1;
  196.     for (i = 1; i <= n1; ++i) {
  197.         if (i < j) {
  198.             mm[i] = j-1;    /* '-1' for shifting to 0-origin    */
  199.             mm[j] = i-1;    /* '-1' for shifting to 0-origin    */
  200.         }
  201.  
  202.         for (k = n / 2; k < j; k /= 2)
  203.             j -= k;
  204.         j += k;
  205.     }
  206.  
  207. /* ---------- Length two butterflies ------------------------------------- */
  208.  
  209.     for (is = 1, id = 4; is < n; is = 2 * id - 1, id = 4 * id) {
  210.         for (i0 = is; i0 <= n; i0 += id) {
  211.             i1 = i0 + 1;
  212.  
  213.             smad (i0, 0);
  214.             smad (i1, 0);
  215.             fft2 (mm[i0], mm[i1]);
  216.             ++ad[i0];
  217.             ++ad[i1];
  218.             cntx[0] += 2;
  219.             cntx[1] += 2;
  220.             cntx[2] += 2;
  221.             cntx[3] += 2;
  222.         }
  223.     }
  224.  
  225. /* ---------- L shaped butterflies --------------------------------------- */
  226.  
  227.     n2 = 2;
  228.     for (k = 2; k <= m; ++k) {
  229.         ++mad;
  230.         n2 = n2 * 2;
  231.         n4 = n2 / 4;
  232.         n8 = n2 / 8;
  233.         e = 6.2831853071719586 / n2;        /* 2 * pi */
  234.         cc1 = (SC15/2) * 0.7071067811865475;    /* sqrt (0.5) */
  235.  
  236.         for (is = 0, id = n2 * 2; is < n; is = 2 * id - n2, id *= 4) {
  237.             for (i = is; i <= n - 1; i += id) {
  238.                 i1 = i + 1;
  239.                 i2 = i1 + n4;
  240.                 i3 = i2 + n4;
  241.                 i4 = i3 + n4;
  242.  
  243.                 smad (i1, 0);
  244.                 smad (i3, 1);
  245.                 smad (i4, 1);
  246.                 fft3 (mm[i1], mm[i3], mm[i4]);
  247.                 ad[i1] += 1;
  248.                 ad[i3] += 2;
  249.                 ad[i4] += 2;
  250.                 cntx[0] += 3;
  251.                 cntx[1] += 3;
  252.                 cntx[2] += 5;
  253.                 cntx[3] += 4;
  254.  
  255.                 if (n4 != 1) {
  256.                     i1 = i1 + n8;
  257.                     i2 = i2 + n8;
  258.                     i3 = i3 + n8;
  259.                     i4 = i4 + n8;
  260.  
  261.                     smad (i1, 1);
  262.                     smad (i2, 0);
  263.                     smad (i3, 1);
  264.                     smad (i4, 1);
  265.                     fft4 (mm[i1], mm[i2], mm[i3], mm[i4],
  266.                         cc1);
  267.                     ad[i1] += 2;
  268.                     ad[i2] += 1;
  269.                     ad[i3] += 2;
  270.                     ad[i4] += 2;
  271.                     cntx[0] += 4;
  272.                     cntx[1] += 4;
  273.                     cntx[2] += 2;
  274.                     cntx[3] += 6;
  275.                     cntx[4] += 2;
  276.                 }
  277.             }
  278.         }
  279.  
  280.         for (j = 2, a = e; j <= n8; ++j, a += e) {
  281.             a3 = 3 * a;
  282.             cc1 = (SC15/2) * cos (a);
  283.             ss1 = (SC15/2) * sin (a);
  284.             cc3 = (SC15/2) * cos (a3);
  285.             ss3 = (SC15/2) * sin (a3);
  286.  
  287.             for (is = 0, id = n2 * 2; is < n;
  288.                         is = 2 * id - n2, id *= 4) {
  289.                 for (i = is; i <= n - 1; i += id) {
  290.                     i1 = i + j;
  291.                     i2 = i1 + n4;
  292.                     i3 = i2 + n4;
  293.                     i4 = i3 + n4;
  294.                     i5 = i  + n4 - j + 2;
  295.                     i6 = i5 + n4;
  296.                     i7 = i6 + n4;
  297.                     i8 = i7 + n4;
  298.  
  299.                     smad (i1, 0);
  300.                     smad (i2, 0);
  301.                     smad (i3, 2);
  302.                     smad (i4, 2);
  303.                     smad (i5, 0);
  304.                     smad (i6, 0);
  305.                     smad (i7, 1);
  306.                     smad (i8, 1);
  307.                     if (ind = (ad[i3] < (mad-1))) {
  308.                         ++ad[i3];
  309.                         ++ad[i4];
  310.                         cntx[2] += 1;
  311.                     }
  312.                     fft5 (mm[i1], mm[i2], mm[i3], mm[i4],
  313.                         mm[i5], mm[i6], mm[i7], mm[i8],
  314.                         ss1-cc1, -ss1-cc1, cc1,
  315.                         ss3-cc3, -ss3-cc3, cc3, ind);
  316.                     ad[i1] += 1;
  317.                     ad[i2] += 1;
  318.                     ad[i3] += 2;
  319.                     ad[i4] += 2;
  320.                     ad[i5] += 1;
  321.                     ad[i6] += 1;
  322.                     ad[i7] += 2;
  323.                     ad[i8] += 2;
  324.                     cntx[0] += 8;
  325.                     cntx[1] += 8;
  326.                     cntx[2] += 4;
  327.                     cntx[3] += 18;
  328.                     cntx[4] += 6;
  329.                 }
  330.             }
  331.         }
  332.     }
  333.  
  334. /* ---------- force final adjustment ------------------------------------- */
  335.  
  336.     ++mad;
  337.     for (i = 1; i <= n; ++i)
  338.         smad (i, 0);
  339.  
  340. /* ---------- Power spectrum (no sqrt) ----------------------------------- */
  341.  
  342.     fft7 (0, mm[1]);
  343.  
  344.     for (i = 1; i < n/2; ++i)
  345.         fft8 (i, mm[i+1], mm[n-i+1]);
  346.  
  347.     fft7 ((n/2), mm[n/2+1]);
  348.  
  349. /* ---------- end and show stats ----------------------------------------- */
  350.  
  351.     end_fft (file_name, ep_name);
  352.  
  353.     for (i = 0; i < 9; ++i)
  354.         fprintf (stderr, " %5u", i);
  355.     fprintf (stderr, "\n  Load Store Shift   Add  Mult   Adj\n");
  356.     for (i = 0; i < 9; ++i)
  357.         fprintf (stderr, " %5u", cntx[i]);
  358.     fprintf (stderr, "\n");
  359.  
  360.     exit (0);
  361. }
  362.